module ghsrnd
implicit none

contains
function vsktdrnd(eta,beta,T,n)
use gamrnds
use ghs
use linear_operators
integer T,n,i,j
double precision eta,beta(n)
double precision vsktdrnd(T,n),c,h(T),a,normb,&
&ups(1:n,1:n),vr(T*n),ji(T)

a=2.0d0*eta/(1.0d0-4.0d0*eta)
normb=dot_product(beta,beta)

ji=gamrnd(.5d0/eta,.5d0,T)

call drnnor(T*n,vr)
forall(i=1:T)
	h(i)=(1.0d0-2.0d0*eta)/(eta*ji(i))
end forall

if (dot_product(beta,beta).gt.1.0d-8) then
	c=(-1.0d0+dsqrt(1.0d0+4.0d0*a*normb))/(2.0d0*a*normb)
	ups=eye(size(beta))+((c-1.0d0)/normb)*(matvec(beta).xt.matvec(beta))
	ups=.t.chol(ups)
else 
	c=1.0d0-a*normb+2.0d0*((a*normb)**2.0d0)-5.0d0*((a*normb)**3.0d0)+14.0d0*((a*normb)**4.0d0)
	ups=eye(size(beta))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(beta) .xt. matvec(beta))
	ups=.t.chol(ups)
end if


forall(i=1:T)
vsktdrnd(i,:)=-c*beta+c*h(i)*beta+dsqrt(h(i))*matmul(ups,vr(n*(i-1)+1:n*i))
end forall
end function vsktdrnd

function stdrnd(eta,n)
!symmetric Student t
use gamrnds
use linear_operators
integer n
double precision eta
double precision stdrnd(n),h(1),a,r(n)

call drnnor(n,r)
if (eta==0.0d0) then
	h(1)=1.0d0
else
	a=2.0d0*eta/(1.0d0-4.0d0*eta)
	h=(1.0d0-2.0d0*eta)/(eta*gamrnd(.5d0/eta,.5d0,1))
end if
stdrnd=dsqrt(h(1))*r
end function stdrnd

function sktdrnd(eta,beta)
use gamrnds
use ghs
use linear_operators
double precision eta,beta(:)
double precision sktdrnd(1:size(beta)),c,h(1),a,normb,&
&ups(1:size(beta),1:size(beta)),r(1:size(beta))

a=2.0d0*eta/(1.0d0-4.0d0*eta)
normb=dot_product(beta,beta)


h=(1.0d0-2.0d0*eta)/(eta*gamrnd(.5d0/eta,.5d0,1))

if (dot_product(beta,beta).gt.1.0d-3) then
	c=(-1.0d0+dsqrt(1.0d0+4.0d0*a*normb))/(2.0d0*a*normb)
	ups=eye(size(beta))+((c-1.0d0)/normb)*(matvec(beta).xt.matvec(beta))
	ups=.t.chol(ups)
else 
	c=1.0d0-a*normb+2.0d0*((a*normb)**2.0d0)-5.0d0*((a*normb)**3.0d0)+14.0d0*((a*normb)**4.0d0)
	ups=eye(size(beta))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(beta) .xt. matvec(beta))
	ups=.t.chol(ups)
end if

call drnnor(size(beta),r)
sktdrnd=-c*beta+c*h(1)*beta+dsqrt(h(1))*matmul(ups,r)
end function sktdrnd
!!
function ghstdrnd(eta,si,beta)
use ghs
use bessels
use linear_operators
double precision eta,si,nu,gam,beta(:)
double precision ghstdrnd(1:size(beta)),c,h,a,normb,&
&ups(1:size(beta),1:size(beta)),r(1:size(beta))

nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si
c=ghcfun0(beta,eta,si)
h=gam/(gigrnd(-nu,gam,1.0d0)*besselr(nu,gam))

if (dot_product(beta,beta).gt.1.0d-3) then
	ups=eye(size(beta))+((c-1)/dot_product(beta,beta))*(matvec(beta).xt.matvec(beta))
	ups=.t.chol(ups)
else 
	a=besseld(nu+1.0d0,gam)-1.0d0
	normb=dot_product(beta,beta)
	ups=eye(size(beta))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(beta) .xt. matvec(beta))
	ups=.t.chol(ups)
end if

call drnnor(size(beta),r)
ghstdrnd=-c*beta+c*h*beta+dsqrt(h)*matmul(ups,r)
end function ghstdrnd
!!
function ghstdrnd_seed(eta,si,beta,gigvar,normvec)
use ghs
use bessels
use linear_operators
double precision eta,si,nu,gam,beta(:),gigvar,normvec(:)
double precision ghstdrnd_seed(1:size(beta)),c,h,a,normb,&
&ups(1:size(beta),1:size(beta))

nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si
c=ghcfun0(beta,eta,si)
h=gam/(gigvar*besselr(nu,gam))

if (dot_product(beta,beta).gt.1.0d-3) then
	ups=eye(size(beta))+((c-1)/dot_product(beta,beta))*(matvec(beta).xt.matvec(beta))
	ups=.t.chol(ups)
else 
	a=besseld(nu+1.0d0,gam)-1.0d0
	normb=dot_product(beta,beta)
	ups=eye(size(beta))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(beta) .xt. matvec(beta))
	ups=.t.chol(ups)
end if

ghstdrnd_seed=-c*beta+c*h*beta+dsqrt(h)*matmul(ups,normvec)
end function ghstdrnd_seed

!!
function gigrnd(nu,d,gam)
double precision gigrnd,nu,d,gam,b,k,h

b=d*gam
k=d/gam
h=nu


if (nu>0) then
        gigrnd=k*gigru(h,b)
elseif (nu<0) then
        gigrnd=k/gigru(abs(h),b)
endif
end function gigrnd

function gigru(h,b)
use dfport
double precision gigru,h,b
!/******************************************************************
! *                                                                *
! *  Generalized Inverse Gaussian Distribution - Ratio of Uniforms *
! *                                                                *
! ******************************************************************
! *                                                                *
! * FUNCTION :   - gigru samples a random number from the          *
! *                reparameterised Generalized Inverse Gaussian    *
! *                distribution with parameters l > 0 and b > 0    *
! *                using Ratio of Uniforms method without shift    *
! *                for  l <= 1  and  b <= 1  and shift at mode  m  *
! *                otherwise.                                      *
! * REFERENCES : - J.S. Dagpunar (1989): An easily implemented     *
! *                generalized inverse Gaussian generator,         *
! *                Commun. Statist. Simul. 18(2), 703-710.         *
! *              - K. Lehner (1989): Erzeugung von Zufallszahlen   *
! *                fuer zwei exotische Verteilungen,               *
! *                Diplomarbeit, Techn. Universitaet Graz.         *
! * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with    *
! *                 unsigned long integer *seed.                   *
! *                                                                *
! * Implemented by K. Lehner, 1990                                 *
! * Revised     by F. Niederl, August 1992                         *
! ******************************************************************/


double precision:: h_setup
double precision:: b_setup
double precision,parameter:: drittel =  0.3333333333333333  !1/3
double precision,parameter:: pdrittel = 0.037037037037037   !1/27
double precision,parameter:: pi=3.14159265358979
double precision m,linvmax,vminus,vdiff,b2,hm12,a,d,e,c;
double precision r,s,t,p,q,eta,fi,fak,y1,y2,max,invy1,invy2,vplus,hm1,xm,ym
double precision u,v,x,z

h_setup=-2.0
b_setup=-2.0

 if ((h.le. 1.0) .and. (b .le. 1.0)) then
				  !/* NO SHIFT m */

	if ((h .ne. h_setup) .or. (b .ne. b_setup)) then
	 
				  !/* Set-up */
	  e = b*b
	  d = h + 1.0
	  ym = (-d + dsqrt(d*d + e))/b
	  d = h - 1.0
	  xm = (d + dsqrt(d*d + e))/b
	  d = 0.5*d
	  e = -0.25*b
	  r = xm + 1.0/xm
	  s = xm*ym
							!/* a = vplus/uplus */
	  a = dexp(-0.5*h*dlog(s) + 0.5*dlog(xm/ym) - e*(r - ym - 1.0/ym))
						!/* c = 1/dlog{dsqrt[hx(xm)]} */
	  c = -d* dlog(xm) - e*r
								!  /* vminus = 0 */

	  h_setup = h
	  b_setup = b
    endif !/* End - Setup */

				  !/* Generator */
	 u = drand1(0) !                                        /* U(0/1) */
	 v = drand1(0) !                                       /* U(0/1) */
	 x = a*(v/u)
                  
    do while (((dlog(u)) > (d*dlog(x) + e*(x + 1.0/x) + c))) !/* Acceptance/Rejection */
	
	 u = drand1(0)                                        !/* U(0/1) */
	 v = drand1(0)                                        !/* U(0/1) */
	 x = a*(v/u)
    end do
else
				  !/* SHIFT BY m */

   	 if ((h .ne. h_setup) .or. (b .ne. b_setup)) then
	 

				  !/* Set-up */
	  hm1 = h - 1.0
	  hm12 =hm1*0.5
	  b2 = b*0.25
	  m = (hm1 + dsqrt(hm1*hm1 + b*b))/b                 !/* Modus      */
	  max = dexp(hm12*dlog(m) - b2*(m + (1.0/m)))        !/* dsqrt[hx(m)] */
	  linvmax = dlog(1.0/max)

					!/* Find the points x1,x2 (-->invy1,invy2) where
					!the hat function touches the density f(x)    */
	  r = (6.0*m + 2.0*h*m - b*m*m + b)/(4.0*m*m)
	  s = (1.0 + h - b*m)/(2.0*m*m)
	  t = b/(-4.0*m*m)
	  p = (3.0*s - r*r)*drittel
	  q = (2.0*r*r*r)*pdrittel - (r*s)*drittel + t
	  eta = dsqrt(-(p*p*p)*pdrittel)
	  fi = dacos(-q/(2.0*eta))
	  fak = 2.0*dexp(dlog(eta)*drittel)
	  y1 = fak*dcos(fi*drittel) - r*drittel
	  y2 = fak*dcos(fi*drittel + 2.0*drittel*pi) - r*drittel
	  invy1 = 1.0/y1
	  invy2 = 1.0/y2

	  vplus = dexp(linvmax + dlog(invy1) + hm12*dlog(invy1 + m)-b2*(invy1 + m + 1.0/(invy1 + m)))
	  vminus = -dexp(linvmax + dlog(-invy2) + hm12*dlog(invy2 + m)-b2*(invy2 + m + 1.0/(invy2 + m)))
	  vdiff = vplus - vminus
									!/* uplus = 1 */

	  h_setup = h
	  b_setup = b
      endif !/* End - Setup */

				  !/* Generator */
		  u = drand1(0)                                  !/* U(0/1)   */
		  v = vminus + drand1(0)*vdiff                   !/* U(v-/v+) */
		  z =v/u
         
	    do while (z < (-m))
		  u = drand1(0)                                  !/* U(0/1)   */
		  v = vminus + drand1(0)*vdiff                   !/* U(v-/v+) */
		  z =v/u
        end do
	  x = z + m
                  
	  do while ((dlog(u) > (linvmax + hm12*dlog(x) - b2*(x + 1.0/x))))  ! /* Acceptance/Rejection */
		  u = drand1(0)                                  !/* U(0/1)   */
		  v = vminus + drand1(0)*vdiff                   !/* U(v-/v+) */
		  z =v/u
         
	     do while (z < (-m))
		  u = drand1(0)                                  !/* U(0/1)   */
		  v = vminus + drand1(0)*vdiff;                   !/* U(v-/v+) */
		  z =v/u
         end do
	  x = z + m;
      end do
 endif
 gigru=x
contains

function drand1(u)
integer u
double precision:: drand1,res(1)
call drnun(1,res)
drand1=res(1)
end function drand1
end function gigru

end module ghsrnd